NIK and RRPlots
The libraries
library(survival)
library(FRESA.CAD)
## Loading required package: Rcpp
## Loading required package: stringr
## Loading required package: miscTools
## Loading required package: Hmisc
##
## Attaching package: 'Hmisc'
## The following objects are masked from 'package:base':
##
## format.pval, units
## Loading required package: pROC
## Type 'citation("pROC")' for a citation.
##
## Attaching package: 'pROC'
## The following objects are masked from 'package:stats':
##
## cov, smooth, var
library(Biobase)
## Loading required package: BiocGenerics
##
## Attaching package: 'BiocGenerics'
## The following object is masked from 'package:pROC':
##
## var
## The following objects are masked from 'package:stats':
##
## IQR, mad, sd, var, xtabs
## The following objects are masked from 'package:base':
##
## anyDuplicated, append, as.data.frame, basename, cbind, colnames,
## dirname, do.call, duplicated, eval, evalq, Filter, Find, get, grep,
## grepl, intersect, is.unsorted, lapply, Map, mapply, match, mget,
## order, paste, pmax, pmax.int, pmin, pmin.int, Position, rank,
## rbind, Reduce, rownames, sapply, setdiff, sort, table, tapply,
## union, unique, unsplit, which.max, which.min
## Welcome to Bioconductor
##
## Vignettes contain introductory material; view with
## 'browseVignettes()'. To cite Bioconductor, see
## 'citation("Biobase")', and for packages 'citation("pkgname")'.
##
## Attaching package: 'Biobase'
## The following object is masked from 'package:Hmisc':
##
## contents
## The following object is masked from 'package:miscTools':
##
## rowMedians
op <- par(no.readonly = TRUE)
pander::panderOptions('digits', 3)
pander::panderOptions('table.split.table', 400)
pander::panderOptions('keep.trailing.zeros',TRUE)
if (!require("BiocManager", quietly = TRUE))
{
install.packages("BiocManager")
BiocManager::install("seventyGeneData")
}
## Bioconductor version '3.15' is out-of-date; the current release version '3.17'
## is available with R version '4.3'; see https://bioconductor.org/install
library(seventyGeneData)
data(vanDeVijver)
class(vanDeVijver)
## [1] "ExpressionSet"
## attr(,"package")
## [1] "Biobase"
Getting the clinical data
pdata <- pData(vanDeVijver)
ROC Plots
table(pdata$Posnodes)
n y 151 144
pander::pander(table(pdata$TTMevent))
pmroc <- plotModels.ROC(cbind(pdata$TTMevent,-pdata$C1used),name="NIK",thr= -0.4) ## Using paper threshold
par(op)
Node positive data
Node Negative data
RRPlot Cox Model
timeinterval <- 5 # Five years
h0 <- sum(pdata$TTMevent & pdata$RFS <= timeinterval)
h0 <- h0/sum((pdata$RFS > timeinterval) | (pdata$TTMevent==1))
mcox <- coxph(Surv(RFS,TTMevent)~C1used,pdata)
pander::pander(summary(mcox)$coefficients)
| C1used |
-1.5 |
0.224 |
0.263 |
-5.69 |
1.3e-08 |
index <- predict(mcox,pdata)
rdata <- cbind(pdata$TTMevent,ppoisGzero(index,h0))
RRAnalysisCox <- RRPlot(rdata,atRate=c(0.90,0.95),
timetoEvent=pdata$RFS,
title="NIK: Breast Cancer",
ysurvlim=c(0.00,1.0),
riskTimeInterval=timeinterval)
par(op)
Expected time to event
toinclude <- rdata[,1]==1
obstiemToEvent <- pdata[,"RFS"]
tmin <- 0.1 + min(obstiemToEvent)
tmax <-max(obstiemToEvent)
sum(toinclude)
[1] 101
timetoEvent <- meanTimeToEvent(rdata[,2],timeinterval)
timetoEvent[timetoEvent > 2.0*tmax] <- tmax
lnotime <- log(obstiemToEvent[toinclude])
lnetime <- log(timetoEvent[toinclude])
lmfit <- lm(lnotime~0+lnetime)
sm <- summary(lmfit)
pander::pander(sm)
| lnetime |
0.518 |
0.0375 |
13.8 |
6.07e-25 |
Fitting linear model: lnotime ~ 0 + lnetime
| 101 |
0.755 |
0.657 |
0.653 |
plot(timetoEvent,obstiemToEvent,
col=1+rdata[,1],
xlab="Expected time",
ylab="Observed time",
main="Expected vs. Observed",
xlim=c(tmin,tmax),
ylim=c(tmin,tmax),
log="xy")
lines(x=c(tmin,tmax),y=lmfit$coefficients*c(tmin,tmax),lty=1,col="blue")
txt <- bquote(paste(R^2 == .(round(sm$r.squared,3))))
text(tmin+0.005*(tmax-tmin),tmax,txt,cex=0.7)
text(tmin+0.015*(tmax-tmin),tmax,sprintf("Slope=%4.3f",sm$coefficients[1]),cex=0.7)
legend("bottomright",legend=c("No Event","Event","Linear fit"),
pch=c(1,1,-1),
col=c(1,2,"blue"),
lty=c(-1,-1,1)
)
MADerror2 <- mean(abs(timetoEvent[toinclude]-obstiemToEvent[toinclude]))
pander::pander(MADerror2)
4.75
RRPlot Cox Adjusted Model
This time we will include Lymph node status from pathology report and
Estrogen receptor alpha expression measurement from microarray
mcox <- coxph(Surv(RFS,TTMevent)~C1used*(ESR1 + Posnodes),pdata)
pander::pander(summary(mcox)$coefficients)
| C1used |
-0.403 |
0.668 |
0.629 |
-0.640 |
0.52186 |
| ESR1 |
0.123 |
1.131 |
0.255 |
0.481 |
0.63079 |
| Posnodesy |
-0.305 |
0.737 |
0.217 |
-1.401 |
0.16112 |
| C1used:ESR1 |
-1.913 |
0.148 |
0.739 |
-2.588 |
0.00966 |
| C1used:Posnodesy |
0.378 |
1.460 |
0.583 |
0.649 |
0.51661 |
index <- predict(mcox,pdata)
rdata <- cbind(pdata$TTMevent,ppoisGzero(index,h0))
RRAnalysisAdCox <- RRPlot(rdata,atRate=c(0.90,0.95),
timetoEvent=pdata$RFS,
title="Adjusted: Breast Cancer",
ysurvlim=c(0.00,1.0),
riskTimeInterval=timeinterval)
par(op)
Expected time to event
timetoEvent <- meanTimeToEvent(rdata[,2],timeinterval)
timetoEvent[timetoEvent > 2.0*tmax] <- tmax
lnotime <- log(obstiemToEvent[toinclude])
lnetime <- log(timetoEvent[toinclude])
lmfit <- lm(lnotime~0+lnetime)
sm <- summary(lmfit)
pander::pander(sm)
| lnetime |
0.516 |
0.0377 |
13.7 |
1.29e-24 |
Fitting linear model: lnotime ~ 0 + lnetime
| 101 |
0.761 |
0.651 |
0.648 |
timetoEvent <- meanTimeToEvent(rdata[,2],timeinterval)
plot(timetoEvent,obstiemToEvent,
col=1+rdata[,1],
xlab="Expected time",
ylab="Observed time",
main="Expected vs. Observed",
xlim=c(tmin,tmax),
ylim=c(tmin,tmax),
log="xy")
lines(x=c(tmin,tmax),y=lmfit$coefficients*c(tmin,tmax),lty=1,col="blue")
txt <- bquote(paste(R^2 == .(round(sm$r.squared,3))))
text(tmin+0.005*(tmax-tmin),tmax,txt,cex=0.7)
text(tmin+0.015*(tmax-tmin),tmax,sprintf("Slope=%4.3f",sm$coefficients[1]),cex=0.7)
legend("bottomright",legend=c("No Event","Event","Linear fit"),
pch=c(1,1,-1),
col=c(1,2,"blue"),
lty=c(-1,-1,1)
)
MADerror2 <-c(MADerror2,mean(abs(timetoEvent[toinclude]-obstiemToEvent[toinclude])))
pander::pander(MADerror2)
4.75 and 5.02
Calibrating the index
calprob <- CoxRiskCalibration(mcox,pdata,"TTMevent","RFS")
pander::pander(c(h0=calprob$h0,
Gain=calprob$hazardGain,
DeltaTime=calprob$timeInterval),
caption="Cox Calibration Parameters")
h0 <- calprob$h0
timeinterval <- calprob$timeInterval;
rdata <- cbind(pdata$TTMevent,calprob$prob)
RRAnalysisCalAdCox <- RRPlot(rdata,atRate=c(0.90,0.95),
timetoEvent=pdata$RFS,
title="Cal. NIK: Breast Cancer",
ysurvlim=c(0.00,1.0),
riskTimeInterval=timeinterval)
par(op)
Expected time to event
timetoEvent <- meanTimeToEvent(rdata[,2],timeinterval)
timetoEvent[timetoEvent > 2.0*tmax] <- tmax
lnetime <- log(timetoEvent[toinclude])
lmfit <- lm(lnotime~0+lnetime)
sm <- summary(lmfit)
sm <- summary(lmfit)
pander::pander(sm)
| lnetime |
0.512 |
0.0374 |
13.7 |
1.29e-24 |
Fitting linear model: lnotime ~ 0 + lnetime
| 101 |
0.761 |
0.651 |
0.648 |
plot(timetoEvent,obstiemToEvent,
col=1+rdata[,1],
xlab="Expected time",
ylab="Observed time",
main="Expected vs. Observed",
xlim=c(tmin,tmax),
ylim=c(tmin,tmax),
log="xy")
lines(x=c(tmin,tmax),y=lmfit$coefficients*c(tmin,tmax),lty=1,col="blue")
txt <- bquote(paste(R^2 == .(round(sm$r.squared,3))))
text(tmin+0.005*(tmax-tmin),tmax,txt,cex=0.7)
text(tmin+0.015*(tmax-tmin),tmax,sprintf("Slope=%4.3f",sm$coefficients[1]),cex=0.7)
legend("bottomright",legend=c("No Event","Event","Linear fit"),
pch=c(1,1,-1),
col=c(1,2,"blue"),
lty=c(-1,-1,1)
)
MADerror2 <-c(MADerror2,mean(abs(timetoEvent[toinclude]-obstiemToEvent[toinclude])))
pander::pander(MADerror2)
4.75, 5.02 and 5.13
Comparing Risks
Comparing concordance Index
## Correlation Index
cindex <- RRAnalysisCI$c.index$cstatCI
## Cox Index
cindex <- rbind(cindex,RRAnalysisCox$c.index$cstatCI)
## Adjusted Cox Index
cindex <- rbind(cindex,RRAnalysisAdCox$c.index$cstatCI)
## Adjusted and Calibrated Cox Index
cindex <- rbind(cindex,RRAnalysisCalAdCox$c.index$cstatCI)
rownames(cindex) <- c("CI","Cox","Adj. Cox","Cal. Adj. Cox")
pander::pander(cindex)
| CI |
0.698 |
0.698 |
0.653 |
0.743 |
| Cox |
0.698 |
0.698 |
0.651 |
0.745 |
| Adj. Cox |
0.707 |
0.706 |
0.658 |
0.748 |
| Cal. Adj. Cox |
0.707 |
0.706 |
0.657 |
0.751 |
Comparing Risk Ratios Index
## Correlation Index
RRratio <- c(RR=RRAnalysisCI$keyPoints$RR[1],
LCI=RRAnalysisCI$keyPoints$RR_LCI[1],
UCI=RRAnalysisCI$keyPoints$RR_UCI[1])
## Cox Index
RRratio <- rbind(RRratio,c(RR=RRAnalysisCox$keyPoints$RR[1],
LCI=RRAnalysisCox$keyPoints$RR_LCI[1],
UCI=RRAnalysisCox$keyPoints$RR_UCI[1]))
## Adjusted Cox Index
RRratio <- rbind(RRratio,c(RR=RRAnalysisAdCox$keyPoints$RR[1],
LCI=RRAnalysisAdCox$keyPoints$RR_LCI[1],
UCI=RRAnalysisAdCox$keyPoints$RR_UCI[1]))
## Adjusted and Calibrated Cox Index
RRratio <- rbind(RRratio,c(RR=RRAnalysisCalAdCox$keyPoints$RR[1],
LCI=RRAnalysisCalAdCox$keyPoints$RR_LCI[1],
UCI=RRAnalysisCalAdCox$keyPoints$RR_UCI[1]))
rownames(RRratio) <- c("CI","Cox","Adj. Cox","Cal. Adj. Cox")
pander::pander(RRratio)
| CI |
3.81 |
2.08 |
6.96 |
| Cox |
3.81 |
2.08 |
6.96 |
| Adj. Cox |
3.42 |
1.93 |
6.07 |
| Cal. Adj. Cox |
3.42 |
1.93 |
6.07 |
Comparing logRank values
## Correlation Index
SurvDif <- c(chisq=RRAnalysisCI$surdif$chisq,pvalue=RRAnalysisCI$surdif$pvalue)
## Cox Index
SurvDif <- rbind(SurvDif,c(chisq=RRAnalysisCox$surdif$chisq,pvalue=RRAnalysisCox$surdif$pvalue))
## Adjusted Cox Index
SurvDif <- rbind(SurvDif,c(chisq=RRAnalysisAdCox$surdif$chisq,pvalue=RRAnalysisAdCox$surdif$pvalue))
## Adjusted and Calibrated Cox Index
SurvDif <- rbind(SurvDif,c(chisq=RRAnalysisCalAdCox$surdif$chisq,pvalue=RRAnalysisCalAdCox$surdif$pvalue))
rownames(SurvDif) <- c("CI","Cox","Adj. Cox","Cal. Adj. Cox")
pander::pander(SurvDif)
| CI |
28.1 |
7.97e-07 |
| Cox |
28.1 |
7.97e-07 |
| Adj. Cox |
28.5 |
6.46e-07 |
| Cal. Adj. Cox |
28.5 |
6.46e-07 |
Comparing Sensitivity
## Correlation Index
sensi <- RRAnalysisCI$ROCAnalysis$sensitivity
## Cox Index
sensi <- rbind(sensi,RRAnalysisCox$ROCAnalysis$sensitivity)
## Adjusted Cox Index
sensi <- rbind(sensi,RRAnalysisAdCox$ROCAnalysis$sensitivity)
## Adjusted and Calibrated Cox Index
sensi <- rbind(sensi,RRAnalysisCalAdCox$ROCAnalysis$sensitivity)
rownames(sensi) <- c("CI","Cox","Adj. Cox","Cal. Adj. Cox")
pander::pander(sensi)
| CI |
0.891 |
0.813 |
0.944 |
| Cox |
0.891 |
0.813 |
0.944 |
| Adj. Cox |
0.891 |
0.813 |
0.944 |
| Cal. Adj. Cox |
0.891 |
0.813 |
0.944 |
Comparing Specificity
## Correlation Index
speci <- RRAnalysisCI$ROCAnalysis$specificity
## Cox Index
speci <- rbind(speci,RRAnalysisCox$ROCAnalysis$specificity)
## Adjusted Cox Index
speci <- rbind(speci,RRAnalysisAdCox$ROCAnalysis$specificity)
## Adjusted and Calibrated Cox Index
speci <- rbind(speci,RRAnalysisCalAdCox$ROCAnalysis$specificity)
rownames(speci) <- c("CI","Cox","Adj. Cox","Cal. Adj. Cox")
pander::pander(speci)
| CI |
0.397 |
0.328 |
0.469 |
| Cox |
0.397 |
0.328 |
0.469 |
| Adj. Cox |
0.397 |
0.328 |
0.469 |
| Cal. Adj. Cox |
0.397 |
0.328 |
0.469 |
Comparing O/E
OERatio <- NULL
## Cox Index
OERatio <- rbind(OERatio,RRAnalysisCox$OERatio$estimate)
## Adjusted Cox Index
OERatio <- rbind(OERatio,RRAnalysisAdCox$OERatio$estimate)
## Adjusted and Calibrated Cox Index
OERatio <- rbind(OERatio,RRAnalysisCalAdCox$OERatio$estimate)
rownames(OERatio) <- c("Cox","Adj. Cox","Cal. Adj. Cox")
pander::pander(OERatio)
| Cox |
0.806 |
0.656 |
0.979 |
0.0285 |
| Adj. Cox |
0.825 |
0.672 |
1.002 |
0.0519 |
| Cal. Adj. Cox |
0.839 |
0.683 |
1.019 |
0.0831 |
Comparing O/Acum
OARatio <- NULL
## Cox Index
OARatio <- rbind(OARatio,RRAnalysisCox$OARatio$estimate)
## Adjusted Cox Index
OARatio <- rbind(OARatio,RRAnalysisAdCox$OARatio$estimate)
## Adjusted and Calibrated Cox Index
OARatio <- rbind(OARatio,RRAnalysisCalAdCox$OARatio$estimate)
rownames(OARatio) <- c("Cox","Adj. Cox","Cal. Adj. Cox")
pander::pander(OARatio)
| Cox |
1.15 |
0.936 |
1.40 |
0.165 |
| Adj. Cox |
1.18 |
0.961 |
1.43 |
0.104 |
| Cal. Adj. Cox |
1.01 |
0.823 |
1.23 |
0.920 |
Comparing NetBenefit
NetBen <- NULL
## Cox Index
NetBen <- rbind(NetBen,RRAnalysisCox$keyPoints$NetBenefit)
## Adjusted Cox Index
NetBen <- rbind(NetBen,RRAnalysisAdCox$keyPoints$NetBenefit)
## Adjusted and Calibrated Cox Index
NetBen <- rbind(NetBen,RRAnalysisCalAdCox$keyPoints$NetBenefit)
colnames(NetBen) <- rownames(RRAnalysisCox$keyPoints)
rownames(NetBen) <- c("Cox","Adj. Cox","Cal. Adj. Cox")
pander::pander(NetBen)
| Cox |
0.223 |
0.232 |
0.222 |
0.239 |
0.252 |
-0.02795 |
| Adj. Cox |
0.230 |
0.249 |
0.202 |
0.252 |
0.274 |
0.00617 |
| Cal. Adj. Cox |
0.190 |
0.211 |
0.162 |
0.215 |
0.239 |
0.01354 |
Compare the ROC AUC
pander::pander(pROC::roc.test(RRAnalysisCI$ROCAnalysis$ROC.analysis$roc.predictor,
RRAnalysisAdCox$ROCAnalysis$ROC.analysis$roc.predictor))
DeLong’s test for two correlated ROC curves:
RRAnalysisCI$ROCAnalysis$ROC.analysis$roc.predictor and
RRAnalysisAdCox$ROCAnalysis$ROC.analysis$roc.predictor
| -0.992 |
0.321 |
two.sided |
0.692 |
0.702 |
